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A  SURVEY  OF  SOME  PROBLEMS  AND  RECENT 
RESULTS  FOR  PARAMETER  ESTIMATION  AND 
OPTIMAL  CONTROL  IN  DELAY  AND  DISTRIBUTED 
PARAMETER  SYSTEMS 

by 

H.  T.  Banks 


ABSTRACT 

We  survey  a  number  of  applications  and  problems  motivating  our  current 
efforts  on  numerical  techniques  for  parameter  estimation  in  and  optimal 
control  of  delay  and  partial  differential  equations.  We  then  outline  two 
different  approaches  for  establishing  theoretical  convergence  results  for 
estimation  algorithms.  An  application  of  modal  techniques  to  the  investiga¬ 
tion  of  transport  in  brain  tissue  is  briefly  explained.  A  sketch  of  a  con¬ 
vergence  theory  for  spline  techniques  for  function  space  parameter  estimation 
problems  is  given.  v 


§1.  Introduction. 


In  this  lecture  we  shall  first  present  a  brief  account  of  several  areas 
of  applications  which  have  motivated  our  recent  efforts,  both  theoretical  and 
numerical,  on  approximation  methods  for  estimation  and  control  of  infinite 
dimensional  systems.  We  then  shall  sketch  the  general  theoretical  ideas  we 
have  employed  to  establish  convergence  results  for  related  iterative  schemes. 
Finally  we  return  to  two  of  the  applications  and  illustrate  the  use  of  these 
ideas  by  explaining  in  more  detail  our  investigations  for  these  problems.  As 
we  shall  make  clear,  our  efforts  on  many  of  the  problems  mentioned  below 
involve  joint  endeavors  with  colleagues  and  students.  In  addition  to  a  well- 
deserved  thank  you  to  Richard  Ambrasino,  James  Crowley,  Patti  Daniel,  Mary 
Garrett,  Karl  Kunisch,  and  Gary  Rosen,  we  would  also  like  to  publicly 
acknowledge  E.  Armstrong  (NASA  Langley  Research  Center),  R.  Ewing  and 
G.  Moeckel  (Mobil  Research  and  Development  Corp.),  P.  Kareiva  (Brown  University), 
J.  P.  Kernevez  (Universite  de  Technologie  de  Compiegne) ,  W.  T.  Kyner  (University 
of  New  Mexico),  and  G.  A.  Rosenberg  (V.  A.  Medical  Center,  lT.  N.  M.  School  of 
Medicine)  for  numerous  stimulating  discussions  and  suggestions  which  have 
substantially  affected  the  investigations  of  our  group  at  Brown  University. 

Our  discussions  here  focus  on  a  general  class  of  systems  including 
nonlinear  delay  systems 

x(t)  =  f (a,t,x(t) ,xt>x(t  -  x1) , . . . ,x(t  -  xv))  +  g(t)  , 

(1)  x(9)  =  <j)(0)  ,  -tv  <  e  <  0  , 

Q  —  (Ct>  L  ^  ,  •  .  •  ,  T^)  , 


nonlinear  distributed  parameter  systems 
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(2) 


3^u  3^u  ,  3u  . .  . 

— 2  =  ql  — 2  +  q2  +  q3U  +  » 

3t  <Jx 

u(0,x)  =  q4^(x)  ,  (0,x)  =  q5^(x) 


u(t,0)  =  g1(t)  ,  u(t,l)  =  g2(t) 


of  hyperbolic  type,  and  parabolic  systems  of  the  form 


3u  M1  3  .  .  .  3u.  .  .  ,,  . 

3F  =  kU)  3^  p  x) 3x}  +  q2U  +  ’ 
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u(0,x)  =  q3d>(x)  , 


'ail  ai2  ai3  ai4‘ 
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fgl(t)] 


U(t,0),  -~<t,0),  U(t,l),  ^(t.Dj^  =  ,  .  y 

[82  j 


A  typical  estimation  problem  consists  of  the  inverse  problem  of  finding  the 
vector  parameter  q  ,  given  observations  {£..}  of  the  state  (or  components 
of  the  state)  corresponding  to  known  inputs  g  or  .  A  typical  control 

problem  (for  fixed  parameter  values  q  )  might  consist  of  minimizing  a  given 
payoff  or  cost  functional  subject  to  (1),  (2),  or  (3),  over  some  admissible 
class  of  control  functions  g  or  g^  . 
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§2.  Motivating  Examples 
I.  The  LN^  Wind  Tunnel 

The  liquid  nitrogen  wind  tunnel  (National  Transonic  Facility)  currently 
being  constructed  at  NASA  Langley  Research  Center  is  a  cryogenic  wind  tunnel 
for  which  the  cost  of  liquid  nitrogen  alone  is  estimated  at  $6.5  x  10^  per 
year  of  operation.  The  tunnel  represents  the  latest  advances  in  technology 
in  that  essentially  independent  control  of  Mach  number  and  Reynolds  number 
(~  temperature)  is  an  anticipated  feature.  Schematically,  the  tunnel  can  be 
represented  as  in  Figure  1. 


FAN 


GASEOUS 

NITROGREN 

VENT 


SECTION 


TEST 


LIQUID 

NITROGEN- 

INJECTION 


SECTION 


Figure  1 

The  basic  physical  model  relating  states  such  as  Reynolds  no.,  pressure, 
and  Mach  no.  to  controls  such  as  LN^  input,  GN^  bleed,  and  fan  operation 
involves  a  formidable  set  of  partial  differential  equations  (Navier-Stokes)  to 
describe  fluid  flow  in  the  tunnel  and  test  chamber.  This  model  has,  not 
surprisingly,  proved  to  be  very  unwieldy  computationally  and  probably  cannot 
be  used  directly  in  design  of  sophisticated  control  laws.  (Both  open  loop  and 
feedback  controllers  are  needed  for  efficient  operation  of  the  tunnel — and 
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this  is  clearly  a  desirable  goal.  Given  the  current  estimates  of  costs  of 
operation,  the  funds  from  only  a  1  or  2%  savings  in  operation  costs  would 
support  a  nontrivial  amount  of  related  research  by  scientists  and  engineers!!) 

In  view  of  the  schematic  in  Figure  1,  it  is  not  surprising  that 
engineers  (e.g.  see  [18])  have  proposed  lumped  parameter  models  (the  variables 
representing  values  of  states  and  controllers  at  various  discrete  locations 
in  the  tunnel)  with  transport  delays  to  represent  flow  times  in  sections  of 
the  tunnel.  A  specific  example  is  the  model  (see  [  1  ])  for  the  Mach  no.  (in 
the  test  chamber)  loop  which  to  first  order  is  controlled  by  the  fan  guide 
vane  angle  setting  (in  the  fan  section) — i.e.  M(t)  ~  GV A(t  -  r)  where  7 
represents  a  transport  time  from  the  fan  section  to  the  test  section. 

In  addition  to  the  design  of  both  open  loop  and  closed  loop  controllers, 
parameter  estimation  techniques  will  be  useful  once  data  from  the  completed 
tunnel  is  available  (current  investigations  involve  use  of  data  from  a  4 
meter  scale  model  of  the  tunnel) . 

II.  Enzyme  Tubular  Reactors 

Column  reactors  in  which  enzyme  mediated  chemical  reactions  take  place 
to  produce  a  desirable  product  (or  products)  from  a  given  substrate  (or 
substrates)  are  of  some  importance  because  of  che  numerous  potential  applica¬ 
tions  in  commercial  production  (e.g.,  purification  of  fruit  juices,  proteolytic 
treatment  of  beer,  synthesis  of  antibiotics  and  steroids)  .  Research  on  the 
operation  of  such  reactors  has  been  carried  out  in  the  laboratory  of  D.  Thomas 
at  Universite  de  Techno logie  de  Compiegne  for  several  years.  Mathematical 
models  for  these  processes  (which  involve  reaction,  diffusion,  and  convective 
transport)  range  from  simple  plug-flow  (PF)  models  to  full-fledged  diffusion- 
convoet ion- react  ion  ( DCR)  models  [16],  [19].  In  an  attempt  to  formulate 
models  with  the  desirable  accuracy  exhibited  by  the  DCR  models  (which  are 


computationally  expensive  and  unwieldy  to  use,  especially  on  small  computers) 
but  which  approach  the  simplicity  of  the  PF  models  (which  in  these  applications 
prove  too  inaccurate  in  their  representation  of  qualitative  phenomena  to  be  of 


practical  use),  J.  P.  Kernevez  and  his  colleagues  have  proposed  lumped  parameter 
models  with  delays.  In  these  models  there  are  several  delays  representing 
convective  transport  and  a  number  of  diffusive  transport  mechanisms.  One 
version  of  such  models,  which  are  nonlinear  due  to  certain  reaction  velocity 
terms,  is  discussed  in  some  detail  by  P.  Daniel  in  [12]  where  additional 
references  may  also  be  found.  To  investigate  the  accuracy  and  potential 
usefulness  of  these  models,  efficient  methods  for  parameter  estimation 
(unknown  parameters  include  several  delays  as  well  as  kinetic  constants)  end 
control  techniques  for  nonlinear  delay  systems  are  essential. 


III.  Gas  and  Oil  Exploration  and  Recovery 

a)  Reservoir  Engineering  Problems 

The  importance  of  inverse  or  parameter  estimation  problems  in  the  gas  and 
oil  industry  is  rather  well-documented.  One  class  of  problems  [8  J,  [IS], 
[26]  involves  use  of  the  flow  equations  in  a  porous  medium  (a  reservoir  or 
oil/gas  field)  to  determine  the  field  porosity  <t>  (the  ratio  of  pore  volume 
to  total  volume)  and  field  permeability  function  k  .  A  greatly  simplified 
model  would  be  based  on  an  equation  (derived  from  conservation  of  mass  and 
Darcy's  law — see  [11],  [20])  for  the  pressure  p  =  p(t,x,v)  in  a  vertically 
homogeneous  field  of  depth  h  ,  say 


$ch 


9p  _3_.hk  3£ 
9t  9x  u  3x 


)  +  f  , 


where  y  =  fluid  viscosity,  c  =  fluid  compressibility,  and  f  is  a  general 
sink/source  term.  The  field  usually  contains  a  n amber  of  wells  (for  production 


L 
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or  observation  or  both)  and  a  typical  problem  is  to  estimate  and  k  or, 

alternatively,  the  total  pore  volume  0  =  i  ^  dh  ,  from  observations  of  p  at 

)  j 

the  well  heads. 

More  realistic  models  involve  several  miscible  fluids  and 
coupled  set  of  partial  differential  equations  [14],  but  the  funca-u.  -itr 1  inverse 
problem  is  similar,  only  much  more  complicated,  of  course, 
b)  Seismic  Exploration 

A  second  class  of  inverse  problems  concerns  determination  of  the  elastic 
properties  of  an  inhomogeneous  medium  via  surface  observations  after  perturbing 
"shocks"  have  produced  waves  in  the  medium.  Models  usually  involve  the 
equations  of  elasticity  [  2  ],  [17];  for  example,  in  the  1-dimensional  problem 
one  might  consider 

?t 

where  p  is  the  mass  density  and  E  =  A  +  2u  for  compressional  or  P-waves, 

E  *  p  for  shear  or  S-waves  with  A,p  the  Lam<5  parameters.  The  boundary 

conditions  at  z  =  0  (here  z  is  the  vertical  distance  from  the  surface) 

include  excitation  or  perturbation  of  the  medium  (often  this  source  input 

itself  is  a  quantity  to  be  "identified") .  From  observations  at  the  surface 

2*0  (these  observations  tisua  1  1  v involve  the  unknown  source  input  and  a 
3u 

velocity  term  for  particle  displacement) ,  one  wishes  to  determine  the 

unknown  functions  p  and  E  and,  in  addition,  the  source  term  if  it  is 

unknown . 

IV.  Large  Space  Structures 

Another  class  of  control  and  identification  problems  for  which  the  models 
are  based  on  the  equations  for  elastic  structures  are  those  dealing  with 


large  space  antennas.  One  such  antenna  that  is  currently  being  developed  bv 
NASA  is  the  Maypole  Hoop/Column  antenna  which  is  depicted  in  Figure  2. 


Figure  2 

This  antenna,  which  when  fully  deployed  somewhat  resembles  an  inverted 
umbrella  (100  m.  in  diameter) ,  consists  of  a  membraneous  surface  of  gold- 
plated  molybdenum  reflective  mesh,  a  collapsible  hoop  or  ring  on  which 
surface  is  stretched,  and  a  telescoping  column  to  which  the  antenna  surface 
is  anchored  and  on  which  feed  assemblies  are  mounted.  The  antenna  in  col¬ 
lapsed  configuration  (similar  to  the  popular  "travel  umbrellas"  that  collapse 
to  fit  into  a  briefcase  or  small  suitcase)  is  to  be  transported  into  space 
in  the  space  shuttle;  it  is  then  deployed  for  use  as  a  communication  antenna. 
The  antenna  surface  itself  is  flexible  and  its  shape  (and  hence  focusing 
properties)  can  be  changed  via  control  stringers  attached  to  48  equally 
spaced  radial  teflon  coated  graphite  "cords"  (  4  control  stringers  per 
radial  cord  or  "gore"  edge).  In  addition  to  the  dynamic  identification  and 
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control  problems  associated  with  initial  deployment  o^  the  Maypole  Hocp/Colunn, 
it  is  anticipated  that  after  long  periods  of  operation,  the  reflector  surface 
will  (due  to  changes  in  elastic  properties  and  forces)  require  adjustment. 

Thus  a  static  problem  of  interest  consists  of  the  following:  Determine  from 
observations  (through  sensing  devices  placed  on  the  gore  edges  on  the  surface) 
the  present  configuration  of  the  antenna  surface  and  then  effect  the  desired 
configuration  or  "displacement"  through  adjustment  of  the  control  stringers. 

A  typical  problem  then  might  involve  a  partial  differential  equation  for 
the  displacement  of  a  circular  membrane  or  thin  plate,  say  P(D,q)u  =  f  , 
where  D  represents  spatial  differential  operators,  q  represents  elastic 
parameters  to  be  estimated,  and  f  entails  applied  forces.  A  simple  example 
might  be 


3  ,  v  3u.  ,  3  ,E  3u. 

3?(rE3^  +  36 (7 


f 


where  E  =  E(r,9)  . 

V.  Dispersion  Models  in  Ecology 

An  important  problem  to  population  ecologists  [21],  [22]  concerns  the 
movement  of  insects  (or,  more  generally,  herbivores)  through  vegetation 
patches.  Outbreaks  or  cyclical  population  explosions  of  some  insects  are 
observed  and  it  is  believed  that  the  nature  of  the  transport  mechanisms  for 
the  insects  affect  the  occurrence  (or  lack  thereof)  of  outbreaks  and  their 
magnitude  and  periodicity.  Typical  equations  to  describe  movement  of  the 
insects  involve  both  diffusive  and  advective  (convective)  terms,  e.g.  for 


I?  *  >>  -  &  « *  '<■> 


1-dimensional  models 
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in  addition  to  the  usual  sink/source  terms  f  .  Depending  on  the  species 
involved,  it  is  generally  expected  that  D  and/or  V  can  depend  on  N'  ,  the 
population  level,  and/or  x  ,  the  spatial  variable.  The  diffusion  coefficient 
D  and  the  convective  velocity  V  may  alternatively,  or,  in  addition,  dene-d  on 
temperature  or  time  (e.g.  as  in  seasonal  migration  of  pests). 

There  are  numerous  estimation  and  control  problems  of  importance  in  the 
context  of  ecological  investigations.  Typically  one  wishes  to  determine  the 
coefficient  functions  D  and  V  from  observations  of  N  and  once  this  is 
done,  one  might  wish  to  estimate  the  optimal  vegetation  density  in  a  patch 
in  order  to  hold  population  levels  in  the  patch  to  a  minimum,  or  at  least 

below  some  given  level. 

» 

VI.  Transport  Models  in  Physiology 

In  physiology  a  great  deal  of  research  is  devoted  to  questions  concerning 
transport  mechanisms  such  as  simple  (passive)  diffusion,  bulk  flow  or  convective 
transport,  facilitated  diffusion,  and  active  transport.  An  example  is  the 
effort  [9  ],  [10],  [23],  [24]  devoted  in  recent  years  to  the  controversy 
involving  bulk  flow  vs.  molecular  diffusion  of  brain  interstitial  fluid  in 
gray  and  white  matter.  The  mathematical  models  again  are  based  on  the 
convection-diffusion  equation 


8u 

3t 


+  V 


8u 

8x 


for  the  concentration  u  of  a  labeled  substance  such  as  sucrose  in  brain 
tissue.  From  experimental  data  (for  u  at  various  times  and  locations  in  the 
tissue  samples)  one  seeks  to  estimate  values  for  V  and  D  in  gray  and  white 
matter  and  contrast  the  transport  properties  of  each  type  of  tissue.  We  shall. 


in  a  subsequent  section  of  this  presentation,  discuss  in  some  detail  an 
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application  to  these  transport  problems  of  some  of  the  methods  that  are  the 
focus  of  attention  in  this  lecture. 


§3.  Theoretical  Foundations 

We  turn  now  to  a  discussion  of  the  theoretical  techniques  that  one  can 
employ  to  establish  convergence  results  for  certain  approximation  schemes  for 
nonlinear  systems  such  as  (1),  (2)  or  (3).  For  the  sake  of  brevity  we  shall 
restrict  our  considerations  to  parameter  estimation  problems.  A  discussion 
of  the  use  of  the  ideas  presented  here  in  control  problems  can  be  found  in 
[12]  in  the  case  of  nonlinear  delay  systems  while  the  cose  of  distributed 
parameter  control  problems  is  considered  briefly  in  [  6  ] . 

For  the  purposes  of  illustration  we  shall  use  a  least  squares  formulation 
(for  a  discussion  of  maximum  likelihood  estimator  ideas,  see  [  3  ]  and  the 
references  therein)  of  the  parameter  estimation  problem.  In  particular,  one 
seeks  to  minimize 

J(q)  =  ~  Z  ! y < t . ;q)  -  C. I2 
i=l 

over  a  given  set.  Q  of  admissible  parameters.  Here  is  an  observation 

for  the  output  t  -*■  y(t;q)  at  t^  with  y(t)  =  Cx(t)  in  the  case  of  (1), 

y(t)  =  col(Cu(t,x..)  , .  .  .  ,Cu(t,x  )  ,Vu  (t  ,x_  )  , . . .  ,Pu  (t  ,x  ))  in  the  case  of  (2), 
i  p  t  1  t  p 

and  y(t)  =  col(Cu(t,x.. ) , . . .  ,Cu(t,x  ))  in  the  case  of  (3),  where  C  and  V 

1  P 

are  matrix  operators  of  appropriate  dimension  and  rank. 

Our  approach  entails  rewriting  (1),  (2),  or  (3)  as  an  abstract  equation 

z(t)  =  A(q) z(t)  +  G(t) 

(4) 

z(0)  =  zQ 

in  an  appropriately  chosen  Hilbert  space  Z  .  The  operator  A  may  be  linear 
or  nonlinear  and  depends  on  the  unknown  parameters  q  .  We  reformulate  the 


estimation  problems  as  ones  of  minimizing 
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J(q)  =  y  E  |r(z(t  ;q>)  -  U2 
i=l 


where  y(t;q)  =  T(z(t;q))  is  an  appropriately  defined  output. 

We  take  a  classical  Ritz-Galerkin  type  approach  to  reducing  these  infinite¬ 
dimensional  state  space  problems  to  a  sequence  of  approximating  finite  dimen¬ 
sional  state  space  problems  that  are  readily  solved  numerically.  For  a  given 

N  ,  N  S 

sequence  Z  of  "subspaces"  of  Z  with  "projections"  P  :  Z  -*■  Z'  ,  we 

minimize 

/(q)  =  y  E  |r(zN(t  ;q))  -  £  i? 

1=1 

N 

over  Q  whcrt*  z  is  the  solution  of  an  approximating  system 


(5) 


zN(t)  -  AN(q)zN(t)  +  PNG(t) 


N 

z  (0)  -  P  zQ  . 


N  N  N 

In  the  methods  discussed  here  we  always  take  A  (q)  =  P  A(q)P  and  obtain 

N 

a  state  convergence  z  (t;q)  ■+■  z(t;q)  .  The  ultimate  goal,  of  course,  is  to 

— N 

insure  convergence  of  some  sequence  {q  }  of  solutions  of  approximating 
estimation  problems  involving  (5)  to  a  solution  q  of  the  problems  involving 
(4).  This  objective  can  be  attained  in  the  cases  of  the  "modal"  and  "spline” 
schemes  we  have  developed  and  tested  numerically  in  [ 4  ] ,  [ 5  ] ,  [ 6  ] ,  [  /  ] , 


[12]. 

To  date  we  have  employed  two  different  theories  to  establish  state  and 
parameter  convergence.  For  distributed  parameter  systems,  both  modal  [6]  and 
spline  [7]  schemes  have  been  investigated  using  an  abstract  semigroup  formulation 

and  Trotter-Kato  type  theorems.  Briefly,  one  establishes  that  the  linear 
N 

operators  A  and  A  (we  suppress  the  q  dependence  here)  satisfv  a 
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uniform  (in  q  and  N  )  dissipativeness  condition  and  generate  C^-scmigroups 
N 

T(t)  and  T  (t)  respectively.  Then  treating  the  nonlinearities 

G(a)  =  F(q,a,z(a))  (  F  is  defined  in  an  appropriate  manner  using  f  from 

(2)  or  (3) )  as  perturbations,  one  considers  in  place  of  (4)  and  (3)  the  implicit 

equations 


(6) 

and 


z(t) 


T(t) 


T(t  -  a)F(q,o,7(^)d'j 


(7)  zN(t)  =TN(t)PNz0+  f  TN(t -0)PNF(q,c,zN(a))d^  . 

> 0 

The  basic  tool  then  is  a  Trotter-Kato  type  result  which,  under  the  conditions: 


(8i) 


|TN(t)|  <  Me<1)t 


for  some  M  and  w  independent  of  N  ; 


(8ii) 


there  exists  a  set  V  C  Dom(A)  ,  V  dense  in  Z  , 
such  that  (Aq-A)P  is  dense  in  Z  for  some  X^  >  0  ; 


(8iii) 


A^z  -  Az  ]  -*■  0 


for 


z  £  V  ; 


guarantees  the  convergence 

(8iv)  T^(t)z  T(t)z  for  z  £  Z  ,  uniformly  in  t  . 

N 

The  convergence  in  (8iv) ,  along  with  (6),  (7)  and  P  I  strongly,  can 

N 

be  used  to  argue  state  convergence  z  (t)  ■+  z(t)  .  This  in  turn  can  be  used  to 
establish  a  desired  parameter  convergence  (for  the  rather  technical  details — 
which  are  nontrivial  when  the  full  dependence  of  the  operators,  projections, 
semigroups  and,  in  some  cases,  the  subspaces,  on  the  unknown  parameters  q  is 


taken  into  account — one  should  consult  [6  ]  and  [7  ])• 


A  somewhat  different  approach  to  spline  methods  for  delay  systems  (’  )  has 
been  taken  in  [4  j,  [5],  [12]  where  the  nonlinearity  f  is  treated  directly 
as  part  of  a  nonlinear  operator  A  =  A(t)  (which  is  now  possibly  time 
dependent) .  In  this  case  one  uses  the  implicit  equations 


z(t)  =  zQ  + 


(A(o)z(o)  +  G(o) ’da 


M  M  M 

z  (t)  =  P  zQ  +  j  {A  (o)z  (o)  +  P  G(o)}do 


in  place  of  (4)  and  (5) .  Under  reasonable  conditions  on  f  one  can  establish 
dissipative  type  inequalities 


<KAN(o)z  -  AN(0)w,  z  -  w»  < 


O)(0)«Z  -  w,  z 


w» 


where  «  ,  »  is  a  specially  defined  inner  product  on  Z  .  With  some 
elementary  analysis  and  use  of  a  Gronwall  inequality,  one  then  obtains  estimates 
for  |zN(t)  -  z(t)|  in  terms  of  integrals  of  (AN(0)  -  A(o)}z(c)  . 

Desired  convergence  results  then  follow  from  convergence  properties  of  A^  . 

Again  the  technical  details  become  quite  involved  when  one  treats  general  non¬ 
linear  delay  systems  with  multiple  unknown  delays.  These  can  be  found  in 
[  5],  [12]. 

N 

We  remark  that  one  need  not  have  Z  a  subspace  of  Z  in  carrying  out 

the  above  theories.  Indeed  in  both  cases  (delay  systems  with  unknown  delays, 

distributed  parameter  systems  with  unknown  coefficients)  outlined  above,  one 

N  N 

finds  that  the  appropriate  P  ,  Z  ,  and  Z  all  depend  on  the  unknown  parameters 

q  (through  the  domain  of  the  function  space  in  the  case  of  unknown  delays  and 

N 

through  the  inner  product  for  Z  and  Z  in  the  case  of  some  distributed 
parameter  examples  as  well  as  the  unknown  delays  problems)  which  of  course 
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vary  as  one  iterates  through  the  sequence  of  approximating  problems  (i.e.,  on 
N  ).  This  feature  results  in  interesting  difficulties  from  both  a  conceptual 
and  computational  viewpoint. 
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§4.  An  Application  to  Transport  in  Brain  Tissue 


We  return  to  the  example  VI  of  §2  involving  the  transport  of  labeled 
sucrose  in  gray  and  white  matter.  A  complete  description  of  the  experimental 
procedures  and  the  questions  being  investigated  can  be  found  in  [24]-  Briefly 
cats  are  anesthetized  and  experiments  of  either  1,  2,  or  4  hours  duration  are 
carried  out.  Labeled  sucrose  is  perfused  into  the  lateral  ventricle.  At  the 
end  of  the  perfusion  period,  the  animals  are  sacrificed  and  their  brains  are 
removed  and  frozen.  Well-stained  areas  of  gray  and  white  matter  perpendicular 
to  the  ventricular  surface  (along  the  x-axis  in  our  notation  below)  are  sec¬ 
tioned  and  analyzed.  This  yields  data  corresponding  to  a  fixed  time  t^  for 
a  maximum  of  4  spatial  locations  x^ ,j  =  1,...,4  ,  in  gray  matter  and  8 

spatial  locations  x.  in  white  matter.  From  this  data  (u(t.,x.)}  for  the 

1  -^-1 

concentration  u  ,  one  wishes  to  compare  transport  in  gray  matter  with  that  in 
white  matter.  The  primary  questions  pertain  to  transport  via  molecular 
diffusion  alone  vs.  transport  via  diffusion  and  bulk  (convective)  flow.  In 
particular,  the  mathematical  problems  reduce  to  those  of  estimating  D,  V,  an 

co  in 


in  +  v  ^  =  D  i-n 
at  +  v  9x  D  .2 


u(t,0)  *=  CQ  . 


In  the  early  experimental  work,  data  for  only  one  time  t^  (1,  2,  or  4  hours) 
and  for  anywhere  from  4  to  8  spatial  locations  x^  were  available.  A 
substantial  concern  is  whether  one  can  develop  accurate  methods  for  estimation 
of  the  three  parameters  in  question  from  such  limited  data. 

We  have  successfully  applied  the  modal  methods  of  Example  4.4  of  [  6  ]  to 
investigate  these  problems.  We  first  summarize  briefly  the  pertinent  ideas 


behind  the  methods.  For  the  purpose  of  illustration,  consider  the  example 


-I'¬ 


ll  =  q,u  +  q„u 
t  1  xx  m2  x 


u(t  ,0)  =  u(t,l)  =  0 


u(0,x)  =  £(x)  , 


which  can  be  reformulated  in  the  form  (4)  in  Z  =  L^CO.l)  by  choosing 
2 

A(q)  =  q^D  +  q^D  (here  D  is  the  differential  operator  in  L^CO.l)  )  with 
2  1 

Dom(A(q))=  H  n  Hq  .  It  can  be  argued  that  A(q)  is  uniformly  maximal 

dissipative  (although  it  is  not  self-adjoint).  For  the  development  of  modal 

approximation  schemes,  general  spectral  results  given  in  [13]  can  be  employed. 

It  can  be  seen  that  A(q)  is  a  relatively  bounded  perturbation  of  a  discrete 

spectral  operator  and  is  itself  a  discrete  spectral  operator  with  a  resolution 

of  identity  and  associated  eigemnanifolds  and  projections.  The  eigenvalues 

2  2  2 

are  found  to  be  X^ (q)  *  -j  IT  q^  -  q^/2q1  with  associated  eigenfunctions 

(q)  *  exp{-q2x/2q^}sin  jttx  .  The  natural  modes  or  eigenfunctions  form  a 

complete  (but  not  orthogonal)  set  Ln  7.  -  1.^(0, 1)  .  However  a  choice  of 

Z  =  span{^ (q)  , . . .  ,4»N(q)  }  ,  while  desirable  theoretically,  is  not  useful  in 

parameter  estimation  algorithms  since  the  basis  elements  are  then  dependent 

upon  the  unknown  parameters  (and  thus  change  with  each  new  estimate  of  the 

q’s  ).  One  can  use  instead  the  near-modal  functions  0^ (x)  =  fl  sin  j“x  and 

take  ZN  ■  span{^, . . .  ,<I>N}  with,  of  course,  A1'  =  PNAPN  where  PN  is  the 

N 

canonical  projection  of  Z  onto  Z 

Convergence  can  be  argued  using  the  Trotter-Kato  formulation  of  (8i)  -  (Siv) 
above.  The  stability  condition  (8i)  follows  immediately  from  the  uniform 

OO 

-  - 

dissipativeness.  Choosing  V  ■  U  Z  (q)  ,  where  q  is  a  limit  of  the  sequence 
-  N  N=1 

of  estimates  q  ,  the  spectral  results  yield  (8ii)  trivially  while  one  must 
work  somewhat  more  to  establish  (8iii)  . 

With  regard  to  implementation,  the  scheme  offers  some  nice  computational 
features  since  the  matrix  realizations  of  the  operators  A‘  (q)  are  given  by 
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r  .2  2  . 

-q1i  tt  1  =  j  . 

[AN(q)]ij  =  *  0  i  ^  j  ,i+1  even, 

2iqA~~]  i  +  j  ,  i  +  j  odd. 

i  -j 

Turning  to  our  investigation  of  these  methods  for  possible  use  in  the 
brain  transport  questions,  we  first  tested  the  methods  with  an  example  for 
which  the  solution  was  "known".  (J.  Crowley  and  M.  Garrett  carried  out  the 
computations  for  this  problem.  J.  Saltzman  supplied  a  "known"  solution 
technique  involving  an  infinite  series  which  was  used  to  generate  "data" 
corresponding  to  fixed  parameter  values  in  the  equation.  This  technique  is 
totally  unrelated  to  the  methods  we  were  testing.)  The  example  used  was 


2 

where  4>(x)  =  a^x  +  a^x  +  is  a  quadratic  satisfying  6(0)  =  1  ,  6(1)  =  0  , 
and  max  <j>  =  <)>(~r)  •  "Data"  were  generated  corresponding  to  true  values 
q*  =  .3  ,  q*  *  1.75  ,  and  q*  =  1.0  .  A  number  of  numerical  trials  with  Jae 
above  described  "modal"  scheme  were  conducted  in  which  the  inverse  problem  foi 
varying  amounts  of  "data"  was  "solved."  We  summarize  briefly  some  of  our 
findings.  In  the  examples  presented  here,  the  notation  I  =  k  ,  J  =  p  in  an 
example  indicates  that  the  data  set  for  this  test  consisted  of  values 
u(t^,Xj)  ,  i  =  l,...,k,  j  *  l,...,p  (i.e.,  k  *  p  "observations"  were  employed 

in  the  inverse  problem) . 


Example  1:  It  was  assumed  that  q^  was  known  and  an  attempt  to  fit  the  data 
was  made  by  searching  for  and  •  Initial  guesses  (for  each  value  of  N 
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jl  o  NO 

tried)  were  q^’  =  1.0  ,  q2*  =  0.0  .  The  "converged"  values  (corresponding 

to  N  *  8  or  16  in  this  and  all  the  examples  presented  here)  were: 

I  =  1,  J  =  3:  q  =  1.82  q9  =  .83 
I  =  2,  J  =  3:  q  =  .2979  q?  =  1.7557  . 


Example  2 :  This  example  was  exactly  the  same  as  Example  1  except  initial 
NO  NO 

guesses  q^*  =  .75  ,  q),’  -1.0  (somewhat  closer  to  the  "true"  values  than 
those  used  in  Ex.  1)  were  used.  For  I  =  1,  J  =  3  ,  the  results  were 
q  *  .3009  ,  q2  =  1.7529  ,  quite  acceptable  in  this  case. 

Example  3:  We  investigated  the  effect  of  using  increasingly  more  spatial 
points  in  our  data  grid  (i.e.,  1=1  with  J  =  4,5,6  ).  For  initial  guesses 
q^’^  =  .8  ,  q^’^  *  .9  ,  increasingly  better  estimates  were  obtained  as  the 
number  of  spatial  grid  points  increased.  We  obtained: 


I  = 

1, 

J 

=  3: 

«1  * 

.6115 

*2  ' 

1.4903 

I  ■= 

1, 

J 

=  4: 

*1  “ 

.3018 

^2  C 

1.7468 

I  = 

1, 

J 

=  5: 

<1  - 

.2978 

52  = 

1.7492 

I  = 

1, 

J 

=  6: 

<1  = 

.2984 

52  = 

1.7493  . 

It  is  clear  that  4  points  in  the  spatial  grid  yields  an  adequate  amount  of 
data  for  this  example. 

Example  4:  In  this  case  we  wished  to  estimate  q^,q0  and  the  boundary  con- 

N  0  NO  NO 

centration  q^  •  Initial  guesses  were  q^’  =  .8,  q„’  =  .9,  q‘  *  =  .5  . 


The  converged  values  were: 


-20- 


I  -  1, 

J  =  6: 

51  = 

.2990 

e9  =  1.6983 

P3  =  1.0356 

1=2, 

J  =  6: 

<1  = 

.2997 

q2  =  1.7469 

q3  =  1.0012 

Example  5:  As  a  final  test  we  modified  the  initial  function  6  used  in  the 

above  examples  to  <K0)  =  1  ,  <f>(x)  =  0  for  x  ^  0  .  This  represents  the  tvpe 

problem  (one  with  a  discontinuity  in  the  boundary-initial  data)  that  one 

encounters  when  using  the  actual  data  collected  in  the  experiments  with  cats 

described  above.  Again  the  results  obtained  were  encouraging.  With  1=2, 

NO  NO 

J  =  6  and  initial  guesses  q*  *  =  .8  ,  q,’  =  .9  ,  converged  values  of 

=  .3019  ,  =  1.7635  were  found  with  a  residual  sum  of  squares  of 

4.3  x  10~3  . 


In  summary,  the  numerical  tests  reveal  that  it  is  probably  unreasonable 

to  expect  to  solve  with  the  "modal"  methods  the  inverse  problems  for  1=1, 

J  =  3  in  most  cases.  However,  problems  with  I  =  2  ,  J  =  3  correspond  to  a 

reasonable  number  of  spatial  observations  for  the  method  in  some  cases.  (By 

changing  labels  during  the  perfusion  period,  Rosenberg,  Kyner,  and  colleagues 

are  now  collecting  data  with  two  time  grid  points.)  For  d3ta  from  white  matter 

(where  J  *  6  is  feasible),  the  methods  should  prove  useful  in  estimating 

D  ,  V  ,  and  in  the  transport  models. 

We  have,  in  fact,  used  the  methods  with  actual  data  sets  (I  =  1,  J  =  S) 

for  white  matter  supplied  by  Kyner  and  his  associates.  The  "modal"  methods 

appear  to  consistently  perform  in  an  acceptable  manner.  Typical  values  obtained 

—6  ? 

in  solving  the  inverse  problems  are  D  =  2.7  *  10  cmVsec.,  V  =  -5.99  ym/min., 
Cq  =  128.5  ,  values  that  are  consistent  with  expectations  based  on  values 
obtained  by  Kyner  and  associates  using  other  techniques. 

We  anticipate  that  extensions  of  these  methods  (or  perhaps  the  spline 


1 
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methods  developed  in  f  7  ])  will  prove  useful  in  future  investigations  of  the 
channeling  structure  in  white  matter  (in  these  problems  the  velocity  coefficient 
V  will  be  spatially  dependent  as  will  also,  in  some  cases,  the  coefficient  of 


diffusion  D  ). 
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§5.  Estimation  of  Function  Space  Parameters 

The  theory  developed  in  [ 6  ]  and  [  7  )  deals  with  estimation  of  parameters 
in  Euclidean  space  sets.  The  framework  is,  however,  general  enough  to  allow 
one  to  treat  problems  in  which  unknown  function  space  parameters  must  be 
estimated.  In  this  section  of  our  presentation  we  shall  give  a  brief  sketch 
of  hew  one  further  develops  such  a  theory.  At  the  sane  time  we  shall  illustrate 
some  of  the  ideas  fundamental  to  spline  methods  as  opposed  to  the  "modal” 
methods  discussed  earlier. 

In  order  to  demonstrate  the  ideas  we  shall  consider  an  equation  of  the 
porous  media  type  (see  §2.IIIa));  that  is,  we  consider 

(9)  qi(x)f? =  h  (q2(x)!^  +  f 

with  homogeneous  boundary  conditions  u(t,0)  =  u(t,l)  =  0  .  In  relating  this 

to  the  porous  media  application  (then  u  *  pressure) ,  one  might  consider  large 

fields  for  which  the  boundary  terms  are  either  constant  or  slowly  varying  in 

time.  In  either  case,  such  nonhomogeneous  boundary  problems  can  be  transformed 

to  problems  with  homogeneous  boundary  conditions  in  a  quite  standard  manner. 

With  certain  smoothness  assumptions  on  ,q2  ’  t*le  °Perator  in  equation  (9} 

can  be  viewed  as  a  standard  Sturm-Liouville  operator  (i.e.  identify 

q^  ~  k  ,  ~  p  in  the  usual  notation  for  the  coefficient  functions — see  (3) 

above  and  p.  40-42  of  [  6  ]).  For  our  discussions  here  we  shall  assume  that 

q  =  (q^.qj)  is  to  be  chosen  from  a  parameter  space  Q  C  l^O.l)  *  L?(0,1) 

2  3 

satisfying  Q  C  {(q^q^  €  H  *  H  jq^  >  0  ,  0  <  in  <  q^^  <  M}  .  (The  smoothness 
hypothesized  will  guarantee  certain  smoothness  properties  for  the  eigenfunctions 
to  be  discussed  momentarily.) 

We  rewrite  (9)  as  an  equation 


(10) 


z(t)  =  A(q)z(t)  +  F(q,f) 
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in  the  state  space  Z  =  X(q)  =  L  (0,1)  where  we  take  as  inner  product 
fl  2 

*  <J>^q,  .  (Here  the  spaces  do  depend  on  the  unknown  parameters,  a 

q  J0  1  — 

complicating  possibility  we  mentioned  earlier.)  The  operators  in  (10)  are 

given  by  F  =  —  f  and  A(q)tJ/  =  —  D(q_D^)  ,  where  Dom(A(q))  =  n  H"  . 
ql  ql  2 

Simple  integration  by  parts  arguments  yield  <A(q)z,z>  <0  so  that  A (c) 

is  uniformly  dissipative  in  X(q)  .  In  fact  A(q)  is  maximal  dissipative  and 

generates  a  C^-semigroup,  and  we  are  thus  in  a  position  to  consider  (6) ,  (7) 

and  the  Trotter-Kato  approach  to  approximation  schemes. 

To  describe  the  spline  methods  we  need  to  recall  the  definition  of  some 

standard  cubic  spline  basis  elements.  For  any  positive  integer  N  we  let 

=  j/N  ,  j  =  +  3  ,  and  let  BN  ,  j  =  -1,...,N  +  1  ,  be  the  cubic 

J  2  ^ 

spline  that  vanishes  outside  ’  ^as  va^ue  ^  anc*  s^-°Pe  ®  at  ^  * 

N  N 

value  1  and  slope  3N  at  t\  ^  ,  and  value  1  and  slope  -3N  at  t^+^  •  (See  [25j . 
p.  73 — and  note  that  our  elements  here  differ  from  those  of  Schultz  only  by  a 
multiplicative  factor  of  24  .) 

N 

For  our  modified  basis  elements  B.  we  take  the  restriction  to  [0,1]  of 

J 

the  following: 


Bo  =  -  4B-i 
Bi  -  -  4i? 


N  ~N 

B”  =  BV  ,  j  -  2, . . . ,N-2  , 

_N  rN 
B  *  B  -  4B 
N-l  N  N-l 

N  ~N 

bn  =  bn  ~  4bn+i  ' 


We  note  that  these  elements  are  in  Dorn (A (q))  . 


N 

We  define  our  approximation  subspaces  X  (q)  C  X(q)  by 

XN(q)  =  span{B^, . . • ,B?!}  and  let  PN(q)  be  the  canonical  projection  of  X(q) 

0  N  N 

onto  X^(q)  ,  i.e.,  PN(q)i|i  *=  E  <i|i,B^>B^  .  Finally  as  usual  we  take 
AN  =  AN(q)  =  PN(q)A(q)PN(q)  .J  0 

Under  an  assumption  that  Q  is  compact  in  ,  one  can  argue  in 

this  case  that  solutions  to  the  estimation  problems  for  (5)  (or  (7))  do  exist. 

1 

We  in  fact  assume  that  Q  is  compact  in  the  C  x  H*  topology  so  that  we 

henceforth  assume  without  loss  of  generality  (possibly  by  taking  a  subsequence) 

that  we  have  a  sequence  {q^}  of  solutions  to  the  estimation  problems  satisfy- 
N  -  1 

ing  q  -*•  q  in  C  x  H  with  q  £  Q  . 

We  briefly  indicate  the  steps  to  verify  (8i)  -  (8iii)  to  insure  convergence 
of  the  semigroups  generated  by  A‘  (q  )  to  the  semigroup  generated  by  A(q)  . 

(As  we  have  noted  before,  this  is  the  fundamental  convergence  result  needed 
for  both  state  and  parameter  convergence.)  The  stability  requirement  (8i) 
follows  from 


N,  .  ...  N.„N.  N.  N.  N.  .  ^  n 

<A  (q  )z,z>N  =  <A(q  )P  (q  )z,  P  (q  )z>N  **  0  * 

q  q 


N 

the  inequality  being  a  result  of  the  uniform  dissipativeness  of  A(q  )  . 

The  operator  A(q)  has,  by  the  usual  spectral  results,  a  CONS  of  eigen¬ 
functions  (q) }  .  In  (8ii)  we  take 

00 

P  -  u  span{H'1(q),...,'P  (q)}  . 

N=1  N 

It  is  then  easily  seen  that  the  conditions  of  (8ii)  obtain  from  the  completeness 
of  the  and  the  relationship  (A^  -  A(q))4^(q)  »  (A^  -  Aj)'F^(q)  . 

We  finally  consider  (8iii)  and  note  that  the  spaces  X(q)  ,  q  £  Q  ,  are 


-25- 


all  equivalent  (recal)  0  <  1.  **  q^  <  M  )  ,  a  fact  which  plays  a  fundamental  role 
in  the  basic  theory  developed  in  [  6  ].  Indeed  we  may,  in  considering  any 
convergence  results,  equivalently  •  se  the  L ^  topology.  Thus,  to  establish 
(8iii) ,  it  suffices  to  argue 

(ID  AN(qV.  -  A(5)Y 

in  L2  .  From  the  smoothness  assumptions  on  Q  (and  hence  q  )  it  is  easily 

seen  that  ¥  €  .  Since  V.  £  Dom(A(q))  we  also  have  4*  €  .  Hence, 

J  3  2  u 

-  4  1 

it  suffices  to  fix  'p  =  V  (q)  in  H  !"!  and  argue  that  (11)  holds  whenever 

N  -  i 

q  -*•  q  in  C  *  H 

Estimates  similar  to  those  we  need  can  be  found  in  Theorem  6.13,  p.  82  of 

[25].  However  we  cannot  use  those  estimates  directly  since  our  projections 
N  N  N  N 

P  (q  )  (onto  X  (q  ))  are  not  the  same  as  the  standard  projections  (of 
Thm.  6.13)  of  L2  onto  S(N)  *  span{B_^ ,B^, . . . ,BN+1}  .  But  using  fundamental 
ideas  similar  to  those  found  in  [25]  (e.g.,  the  Schmidt  inequality  and  estimates 
for  the  appropriate  interpolating  splines)  one  can  establish: 

For  1)1  6  H4  O  hJ  , 

I*  -  p*Vi  |dSi 

N 

|D(i|>  -  p\)|  <  -f  !d%| 

N 

j D2 (^  -  P%)|  <  -§  |D4<j,I 
N 

N  N  N 

where  the  norms  are  the  usual  norm  ,  P  =  P  (q  )  as  defined  earlier,  and 

the  constants  K^»K2’K3  are  independent  of  N  and  \p  . 

We  thus  have  D^N  -*■  D2^  in  -*■  Di|>  in  C  ,  where  0N  =  PN^  . 

Furthermore  we  know  that  q^  ■*  q^  in  C  ,  q^  q2  in  .  Since  PN  -*•  I 
it  thus  follows  from  elementary  arguments  that 


A(5)i|>  =  (l/q^Dq^ty  +  (q^q^D2^  . 


converges  in  to 

The  assumptions  on  Q  made  above  are  not  unreasonable  for  many 
applications.  We  have  successfully  used  these  spline  methods  in  computational 
packages  for  function  space  parameter  estimation  in  models  for  insect 
dispersion  (see  §2.V).  In  those  applications,  the  smoothness  and  compactness 
assumptions  listed  above  are  satisfied  when  one  formulates  the  problems  and 
parameterizes  Q  in  a  way  that  is  natural  for  and  consistent  with  the 
experimental  and  theoretical  efforts  of  population  ecologists. 
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